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l.O  INTRODUCTION 

The  T6  missile  is  a small  semi-active  terminal  homing  air-to-ground 
missile.  It  consists  of  many  subsystems  such  as  missile  airframe  dynamics, 
fin  actuator  dynamics,  attitude  sensors,  geometric  models,  laser  seeker 
dynamics,  stabilization  filters,  and  guidance  filters.  The  mathematical 
description  of  this  complex  system  results  In  time-varying,  coupled  high 
order,  nonlinear  differential  equations.  A six  degree-of-f reedom  pro- 
gram has  been  written  for  the  simulation  of  this  system.  The  system 
is  currently  modeled  by  high  order,  coupled  linear  scalar  differential 
equations  or  high  degree  transfer  functions  in  the  s domain.  The  stabi- 
lization filter  has  been  designed  to  smooth  the  response  of  the  missile 
and  the  guidance  filter  has  been  designed  to  guide  the  missile  for 
accurate  impact  of  a target. 

Several  performance  studies  of  the  missile  system  have  been  made 

recently.  For  example,  the  continuous-time  linear  model  of  the  missile 

1 2 

system  was  converted  to  a discrete-time  model  ' so  that  a digital 
stabilization  filter  and  a digital  guidance  filter  could  be  synthesized. 
Another  example  is  an  improved  trajectory  which  was  designed  by  using 

3 

modern  dual-model  control  law  theory  so  that  the  missile  Impacts  a tar- 
get at  a predetermined  pitch  attitude  approach  angle. 

In  reviewing  previous  studies,  several  research  problems  have  been 
identified  by  the  authors  for  improving  these  prior  results.  In  this 
report  we  will  concentrate  on  two  subjects:  (I)  accurate  representation 
of  continuous-time  state  equations  of  T6  missile  system  by  discrete-time 
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state  equations  and  (2)  stability  of  coupied  mui tivariable  missile  sys- 
tems. These  two  subjects  were  chosen  due  to  the  possibiiity  for  immediate, 
practical  implementation  and  theoretical  analysis.  For  example,  the 
missile  system  is  currentiy  represented  by  high  order  continuous  state 
equations.  These  state  equations  need  to  be  converted  to  discrete-time 
state  equations  from  which  the  impiementations  of  the  controllers  for 
the  missile  system  can  be  accomplished  by  microprocessors.  In  the  prior 
results,  a crude  approximate  discrete-time  model  has  been  used  for  the 
analysis  and  synthesis  of  digital  autopilot  systems.  Accurate  simulation 
and  control  of  the  missile  system  require  that  a more  accurate  model  be 
used.  In  this  report  we  give  an  approach,  which  is  based  on  a model 
reduction  technique  developed  by  the  authors,  for  the  formulation  of 
accurate  discrete-time  models. 

The  T6  missile  yaw  and  roll  dynamics  are  heavily  coupled  multivariable 
subsystems.  A primary  and  often  vexing  concern  in  the  design  of  coupled 
multivariable  systems  is  the  determination  of  the  stability  of  the  coupled 
systems.  Since  the  methods  for  determining  the  stability  of  coupled  multi- 
variable  systems  are  quite  different  than  those  developed  for  uncoupled 
systems,  an  approach  for  determining  the  stability  of  the  coupled  sys- 
tems has  been  developed  and  is  presented  in  this  report. 
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2.0  ACCURATE  REPRESENTATION  OF  CONTINUOUS-TIME  STATE  EQUATIONS  OF  T6 
MISSILE  SYSTEM  BY  DISCRETE-TIME  STATE  EQUATIONS 

2.1  INTRODUCTION 

The  accurate  description  of  many  practical  systems  results  In  high 
order  continuous-time  state  equations  which  are  difficult  to  work 
with  — whether  one  Is  simulating,  realizing  or  designing  these  systems. 

In  particular,  the  systems  which  consist  of  time-varying  and/or  time- 
delay  components  are  especially  difficult.  However,  the  difficulty  of 
these  complex  processes  can  often  be  significantly  reduced  If  a set  of 
continuous-time  state  equations  can  be  accurately  represented  by  a 
discrete-time  state  equation  set.  This  simplification  can  be  accomplished 
by  several  methods. 

A procedure  often  used  for  small,  time-invariant  systems  Is  to 

directly  evaluate  the  state  transition  matrix  and  its  convolution  with 

the  Input.  This  yields  an  exact  representation  of  a discrete-time  state 

equation.  However,  for  either  a large  system,  a time-varying  system,  or 

a system  with  slight  variations  of  system  parameters,  this  approach  is 

4 5 

not  practical.  Tustin  and  Boxer  and  Thaler  have  Indirectly  constructed 

approximate  discrete-time  state  equations  via  s-  or  z-transforms.  In 

practical  applications,  it  Is  often  necessary  to  determine  directly  an 

approximate  model  in  the  time  domain.  Recently,  Shieh,  et  al^  suggested 

a time  domain  method  for  this  problem  by  using  the  newly  developed  block- 

pulse  functions  (developed  from  Walsh  functions^).  In  this  report,  it 

is  shown  that  the  models  obtained  by  the  referenced  authors^'^'^  are 

equivalent  and  that  this  model  is  one  possible  model  resulting  from  the 

2 

use  of  the  method  presented  here.  A convenient  method  often  used  in 


industry  is  to  construct  an  approximate  discrete-time  state  equation  by 
truncating  the  infinite  series  obtained  from  expansion  of  the  exact 
state  transition  matrix  and  its  convolution  with  the  input  in  the  time 
domain.  However,  the  truncating  error  depends  heavily  on  the  number  of 
terms  and  the  sampling  period  used.  Based  on  the  idea  of  mul ti -feedback 
and  multi -feedforward  control  theory,  a method  for  determination  of  the 
approximate  discrete-time  state  equations  is  presented  in  this  report. 

O 

Matrix  continued  fractions  are  used  as  a basis  for  this  investigation. 

2.2  APPROXIMATION  OF  STATE  TRANSITION  MATRIX 

Consider  the  system  governed  by  the  continuous-time  state  equations, 

X (t)  = Ax  (t)  + Bu  (t) 
o o o 

x^(0)  = x„  (1) 

o o 

where  A and  B are  n X n and  n X p constant  matrices,  respectively,  x^(t) 
is  an  n X 1 column  state  vector,  u^(t)  is  a p X i continuous- time  input 
vector,  and  x^(0)  is  an  initial  vector.  The  solution  of  Eq.  (1)  is 


x^(t)  * 4>(t)  x^(0)  +J  4>(t-X)  Bu^(X)  dX 


where  $(t)  is  the  continuous-time  state  transition  matrix,  or 

$(t)  = e^^  ® [e^^]  = IA'CT)]*^  {2h 

and  t = kT,  k is  an  integer.  If  $(t)  and  the  convolution  Integral  in 
Eq.  (2a)  can  be  explicitly  evaluated,  the  corresponding  discrete-time 
state  equation  is  obtained  by  substituting  t = kT  into  Eq.  (2a).  From 
the  point  of  view  of  practical  implementation,  we  are  interested  in  the 


fW 


|ilfW5P!ljU4  JPM! 


system  in  which  the  input  signal  5s  a piecewlse-constant  and  the 

sampling  instants  are  0,  T,  2T,  which  gives  us  a rectangular  approx- 

imation of  u^(t).  This  is  equivalent  to  inserting  :t  sample-and-zero- 
order-hold  device  before  an  integrator.  The  approximate  input  is 
designated  as  u(t),  or 

u(t)  = u(kT)  a UQ(t)  for  kT  < t < (k+l)T  (3) 

The  approximate  state  due  to  u(t)  is  designated  as  x(t).  The  new  state 
equations  are 

x(t)  = Ax(t)  + Bu(t) 


x(0)  = x^ 


(5*3) 


The  solution  of  Eq.  {ka)  is 


AT 


x(kT  + T)  = e ' x(kT)  + 


£ 


e^**  B da  • u(kT) 


(4b) 


where  a = T - X.  Letting  x(kT)  = x(k)  and  u(kT)  = u(k),  we  have  the 
discrete-time  state  equation 

x(k  + 1)  = $(T)  x(k)  + Lu(k)  (4c) 

where 

$(T)  * e^^ 
and 

L = e''“  B da 

Jo 
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The  solution  of  Eq.  (4c)  is 


u •‘■I 

x(k)  = $(T)  X + X 'l'(k-j-l)  Lu(j) 
j=0 


(4d) 


The  system  matrix  ^(T)  and  its  convolution  with  the  input  can  be 
expressed  by  the  infinite  series: 

$(T)  = = I + AT  + ^(AT)^  + ^(AT)^  + ^(AT)'*  + ... 


= i jr(AT)J 

j=0  J* 


and 


■r 


e^“  B da  = T[l  + ^ AT  + ^(AT)^  + ^(AT)^  + ...]B 


T I TjV  “ 


(5a) 


(5b) 


where  I is  the  identity  matrix.  Eq.  (5b)  can  be  further  simplified  as 


L = T[l  + AT  + |p(AT)^  + ...  - IKAT)"’  B 


[e^^  - I]a"*  B = [$(T)  - I]a"’  B 


(5c) 


For  practical  reasons,  we  are  interested  in  the  approximate  '1>(T), 
denoted  as  ♦*(T)  or  G.  Eq.  (4c)  becomes 


x*(k  + 1)  = G x*(k)  + Mu(k) 


(6a) 
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G = $*(T)  = $(7) 


M = [G  - I]a“’  B a L 

x*(k)  a x(k)  (6b) 

The  solution  of  Eq.  (6a)  Is 

k-1 

x*(k)  = 4'*(k)  x*(0)  + X l-^dt-j-I)  Mu(j)  (6c) 

j=0 

where 

x*(0)  = x{0) 

4>*(k)  = «*(kT)  * G*^ 


k 

where  G is  the  approximate  discrete-time  state  transition  matrix. 

The  natural  question  arises  — how  to  accurately  determine  G and 

2 

M?  A popular  method  often  used  in  industry  is  to  approximate  $(T)  in 
Eq.  (5a)  by  truncating  the  Infinite  series: 


$(T)  = I + AT  + ^(AT)^  + 3y(AT)^  + iy(AT)^  + ...  (7a) 

a I + AT  (7b) 

a I + AT  + j^(AT)^  (7c) 

a I + AT  + J^(AT)^  + Jj-(AT)^  (7d) 

a I + AT  + ~(AT)^  + yj-(AT)^  + |j■(AT)^  (7e) 


If  a sufficiently  large  number  of  terms  In  Eq.  (7a)  Is  used,  then  a 
satisfactory  approximation  may  be  obtained.  In  this  report,  we  introduce 
a geometric  series  to  approximate  the  infinite  series  In  Eq.  (7a)  as 


f o 1 1 ows : 


*(T)  = I + AT  + ^(AT)^  + ^(AT)^  + |j-(AT)^  + ^(AT)^  + ...  (8a) 

CO 

= I + AT  + ^(AT)^  + 2 V(AT)-'  (8b) 

^ j=3-’- 

a I + AT  + y(AT)^  + •^-(AT)^  + -^(AT)^  + ■^■(AT)^  + ...  (8c) 


‘ 2 2'  2 

00 

- I + AT  + ^(AT)^  + Y.  (8d) 

j=3 

= I + [I  + y AT  + •^-(AT)^  + ^(AT)^  + . ..]AT  (8e) 

2 2-’ 

“ I + [I  - j at]"’  at  (8f) 

“ [I  - j at]"’ [I  + i AT]  (8g) 

= $*(T)  (8h) 


The  Infinite  series  In  the  brackets  of  Eq.  (8e)  Is  the  well-known 
geometric  series.  Comparing  Equations  (7c)  and  (8d)  with  Eq.  (8b), 
we  observe  that  the  first  three  terms  In  all  three  equations  are  Iden- 
tical, while  other  terms  differ  by  their  weighting  factors:  zero  In 
Eq-  (7c),  ' In  Eq.  (8d)  and  j(j-l)rj"2)  ..TV  * 

The  approximation  of  4>(T)  In  Eq.  (8a)  given  by  4>*(T)  In  Eq.  (8h) 

Is  much  better  than  that  of  Eq.  (7c). 

For  elaboration  of  the  above  conclusion.  It  Is  convenient  to 
Investigate  the  Inverse  system  of  $*(7),  denoted  by  **(-T).  The 
formulation  of  **(-T)  Is 
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00  ^ ^ 

$*(-T)  = [I  + AT  + y(AT)^  + 2 -!Tr(AT)j’l 

j»3  2-’ 


=[!+[!-  j AT]"’  AT] 

= [I  + [(AT)"’  + (-  i|)]"’]  ’ 

- [H,  + [H2  + [H3]"’]"’]"’  (9a) 

where  = I,  H2  = (AT)  ’ and  ■ -21. 

g 

Equation  (9a)  is  the  formulation  of  a matrix  continued  fraction 
which  is  considered  as  a generalized  mult  I -feedback,  multi -feedforward 
control  system.  A representative  control  system  for  Eq.  (9a)  is  shown 
in  Figure  1.  Note  that  Hj  Is  a feedback  matrix  gain,  while  and  H^’ 
are  the  parallel  feedforward  matrix  gains.  The  use  of  is  very 
important,  for  example,  if  is  neglected  In  Eq.  (9a),  we  have 

$M-T)  = [H,  + [H2]"’]'’ 

- [I  + AT]"’  (9b) 

or 

4>*(T)  = I + AT  (9c) 

Equation  (9c)  is  equal  to  Eq.  (7b),  which  has  the  first  two  terms  of 

the  infinite  series  in  Eq.  (7a).  Therefore  not  only  contributes 
1 2 

the  third  term,  j(AT)  , but  also  generates  an  Important  series, 

00 

j -4— r(AT)-’,  as  shown  in  Eq.  (9a).  From  this  we  conclude  that  If  we 
j“3  2-'"' 

can  establish  more  terms  in  the  inner  loop  of  the  control  system  In 
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Figure  2-1,  we  may  have  a more  accurate  approximation.  For  example,  if 
we  Insert  a constant  matrix  = -3(AT)  ' in  Eq.  (9a),  we  obtain 
«*(-T)  s [H,  + + [Hj  + 

- + (H2  + H^^)1[H,H2H^H^  + (HjH^  + + I]"' 

- [6(AT)‘^  - 2(AT)'’][6(AT)'^  + ^f(AT)"’  + l]"’ 

- [I  - y AT][I  + j AT  + ^ (AT)^]"’ 

- {[I  - y at]"’  [I  + yAT  + y(AT)^]}"’ 

- (I  * at  * ^(at)2  * ^at)3  t ^ T " 

+ ■^(AT)^  + ^^(AT)^  + ...](AT)^  I 
3 y 

- {1  + AT  + ^(AT)^  + ■"  TU^Vi\J  ■ y ^Tl"’(AT)'*r’  (9d) 

or 

$*(T)  = [I  - y AT]"’  [I  + y AT  + y(AT)^] 

= I + AT  + ^ (AT)^  + ^ (AT)^  + Q yj  [I  - 1 AT]"’(AT)'*  (9e) 

A block  diagram  corresponding  to  Eq.  (9d)  Is  shown  In  Figure  2-2.  Hj 

and  Hj  are  mult  I -feedback  gains;  H2  and  are  mult  I -feedforward  gains. 

1 3 

contributes  one  more  Important  term,  yj-(AT)  , than  that  of  in 
Eq.  (9a)  and  another  geometric  series,  (I  - y AT)  ’,  with  the  weighting 
factor  (]  ' Comparing  Eq.  (9e)  with  Eq.  (7d),  It  can  be  seen 

that  Eq.  (9e)  gives  a better  approximation.  These  Hj  . . . can  be 

Q 

generated  by  performing  the  matrix  continued  fraction  expansion  on 
♦(T)  given  In  Eq.  (8a). 
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2.3  APPROXIMATION  BY  MATRIX  CONTINUED  FRACTIONS 

The  inverse  of  the  system  matrix  e^^  in  Eq.  (^c)  can  be  expanded 
into  a matrix  continued  fraction  as  follows: 

e"^^  = i [i  + AT  + |y(AT)^  + 3y(AT)^  + (iOa) 

= IA2J  + A22  + A2J  + ...][Ajj  Aj2  + Ajj  + Ajjij  + ...]  ' (10b) 

= [H,  + [H2  + [Hj  + [H^  + [Hg  + (iOc) 

where 

A21  = 1,  A2^.  - 0 for  j = 2,3.... 

A,,  = I.  Ai  j - .^J^(AT)-'"’  for  j - 2,3.... 

The  matrix  quotients  Hj  in  Eq. (lOc)  can  be  obtained  by  the  matrix  Routh 

q 

array  and  matrix  Routh  algorithm  as  follows: 


'll 


'12 


H,  .A„A-|  < 


A13  A^j^  .., 


"21 


22 


23 


”2  “ ^^21^1  ^ 


”3  “ ^31^'tl  ^ 


^31  “ ^12  “ ”1^22  ^32  " ^13  ' ”/23  ^33 


(Ha) 


\l  “ ^22  ■ ”2^32  \2  “ ^23  " ”2^3 


The  block  elements  of  the  first  and  second  row  of  Eq.  (11a)  are  the 
matrices  In  Eq.  (10b).  The  block  elements  of  the  subsequent  rows  are 
evaluated  by  the  following  matrix  Routh  algorithm: 


f 

1 

) 

l: 

i 


i' 

{ 

\ 

4 


i 
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= 21,  = j(AT)'\  H.^2  “ ^+3  * 2)  (AT)"’ 

for  j = 5,  9.  13,  17 (A  + 4),  ...  (lie) 

The  approximation  of  Eq.  (10)  can  be  obtained  by  truncating  the  matrices 
In  the  inner  positions  of  Eq.  (10c);  i.e., 

= [H,  + + [H3  +[H5  + 

* [H,  + 

* [H,  + + [Hj]"’]"’]"’ 

* [H,  + [H2  + 

* ...  (12) 

By  successively  determining  the  Inversions  of  the  matrices  Hj  In 

Eq.  (12),  we  can  obtain  the  approximation  in  Eq.  (12).  In  order  to 
avoid  these  tedious  inversion  operations,  the  following  technique’®  may 
be  used  to  simplify  the  matrix  inversions. 

To  demonstrate  the  procedure,  we  choose  n ■ 5 with  Hj , j » 1 ...  n 
given.  First,  we  formulate  a chain  of  n matrices. 
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^ - n 


(13a) 


(13b) 


The  structure  of  P Is  as  follows:  The  matrix  quotients  Hj  are  placed 
on  the  diagonal  In  ascending  order;  the  first  matrix  P^  starts  with 

and  the  next,  ?2  with  H2,  etc.  Then  Identity  matrices  are  placed 

on  the  location  above  the  Hj  and  on  the  diagonal  as  shown  In  Eq.  (13b). 

All  other  block  elements  are  zero.  The  product  of  the  chain  of  the 

matrices  In  Eq.  (13a)  Is 

(H,H2HjH^^H5)  (H^H^Hj  + 

0 (H2H3H^Hj) 


(H,  + H3  + Hj)  0 

(H2H3  + H2HJ  + H|,Hj)  I 

(HjH^Hj)  ( 


(Hj  + Hj)  0 

(H^Hj)  I 


-AT 

The  approximation  of  e can  be  obtained  from  the  block  elements  In 
the  first  two  rows  of  P,  or 

e'^'^  a [Hj  + [H2  + [H3  + [H^  + [Hg]"’]'’]'’]"’]'’ 

= + (H^Hj  + H2H5  + + I]- 

+ (H,H2H2  + H,H2H5  + + (H,  + + H^)] 

-II  p, = [i  P2  j p;l][i  p,  j p;l]-’ 

j=i  j=i  '•J  j«i  " j=i  ‘‘j  " 

* ^2,j  ’’ll^ 

j»l  ’ j=l  ’ 


-I 


(14) 


where  P,  . and  P-  . are  the  block  elements  In  the  first  and  second  rows 
of  the  matrix  P,  respectively. 

Using  the  above  procedure  to  determine  P,  , and  P,  and  substl- 
tuting  the  values  of  Hj , yields 


AT 

e 


H,  +[H2  + [Hj  + [H^  + 


n=<x) 


] 


a I + AT 

a [I  - J at]"’ [I  + J AT] 

a [I  - y AT]"’ [I  + y AT  + y(AT)^] 

a [I  - y AT  + -lyCAT)^]  ’[1  + y AT  + •|y(AT)^] 

s [I  - I at  + ^(AT)^]"’[|  + l-AT  + |5-{AT)^  + ^(AT)^] 


(1 5a) 

(15b) 

(15c) 

(15d) 

(I5e) 

(I5f) 

(I5g) 
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Equations  (15d)  and  (15e)  are  equal  to  Eq.  (8g)  and  Eq.  (9e) , respec- 
tively. 

2.4  DERIVATION  OF  DISCRETE-TIME  STATE  EQUATIONS 


Rewriting  Eq.  (6a),  we  have 

x*(k  + i)  * 6x*(k)  + Mu(k)  (i6a) 

where 

G » $*(T)  a «(T)  (i6b) 

M - [G  - I]a"’  B a [$(T)  - |]A~’  B (i6c) 

u(k)  “ u(kT)  a u (kT)  (i6d) 

o 

x*(k)  a x(k)  a Xp(kT)  (l6e) 

The  solution  of  Eq.  (I6a)  is 

x*(k)  = x*(0)  + Y Hu(J)  (17) 

j-0 

The  required  system  matrix  G obtained  In  Eq.  (15)  is 


16 


17 


The  subscript  of  G In  Eq.  (I8)  Indicates  the  number  of  Hj  used.  Substi- 
tuting G obtained  In  Eq.  (I6c)  yields 


M = [G  - I]a"'  B 
5 TB  ^ 

s T[l  - i AT]'’  B = 

3 T[l  - jAT]'’[l  + ^ AT]B  ^ 

= T[l  - i AT  + i^(AT)^]'’B  ^ 

s Til  - |aT  + ^(AT)^]'’ll  + lo  = ”6 

s Til  - IaT  + ^(AT)^  - .^(AT)^]'’ri  + ^(AT)^]B  ^ 

s ...  (19) 

The  subscript  of  M Indicates  the  number  of  Hj  used. 

If  the  polygonal  hold  (which  Is  a device  to  Integrate  by  the 
trapezoidal  approximation”)  can  be  realized,  a more  accurate  discrete- 
time  state  equation  can  be  obtained.  The  staircase  Input  u^(k)  is 

„*(k)  . ,20) 

The  approximate  discrete-time  state  equation  is 

x*(k  + 1)  - G x*(k)  + Mu*(k)  (21) 

where  G and  M are  shown  In  Equations  (18)  and  (19),  respectively. 


Since  more  accurate  inputs  are  used  in  Eq.  (21),  the  response 

x*(k)  in  Eq.  (21)  is  more  accurate  than  the  x*(k)  in  Eq.  (16).  If 

AT 

j = 3 and  Hj,  and  are  used  to  approximate  e , then  and 
should  be  used.  Equation  (21)  then  becomes 

x*(k  + 1)  = [I  - y AT]"’ [I  + i AT]x*(k)  + T[l  - j AT]"’  Bu*(k)  (22) 

The  model  given  in  Eq.  (22)  has  been  derived  by  others^  Therefore, 
it  is  seen  that  their  result  is  a special  case  of  the  results  proposed 
in  this  paper. 

In  the  z-  or  s-domain,  Eq.  (22)  can  be  derived  by  using  the 
transfer  function  of  an  approximate  numerical  differentiator  or 
integrator  : 


Z[l.]  aiilL 


Z[/dt]  a 1 11^ 


■1 


1 


1-z 


1 


(23a) 

(23b) 


In  the  frequency  domain,  the  numerical  integration  operator 
in  Eq.  (23b)  is 

-1 


T(ja))  = j -~ 
^ 1-z 


. T 


u)T 


= -j  J cot  Y- 


(2ka) 


z=e 


jwt 


The  amplitude  and  phase  characteristics  are 


|T(ja))  1 = y |cot  ^1 


/T(ia)) 


TT 

2 


(24b) 


(24c) 
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From  Eq.  (24c)  we  observe  that  we  have  perfect  phase  characteristics 
but  the  amplitude  characteristics  differ  from  the  ideal  characteristics  at 
higher  frequencies.  Note  that  the  system  matrix  in  Eq.  (22)  is  con- 

structed by  using  an  odd  number  of  Hj  (j  = i,2,3).  If  an  odd  number  of 
are  used,  the  approximate  models  of  Gj  give  ideal  phase  characteristics  and 
better  amplitude  characteristics.  For  example,  the  z-transform  of 
Equations  (4a)  and  (21)  with  G = G^  and  M = gives 


Z[x(kT)]  = AZ[x(kT)]  + BZ[u(kT)]  = Ax(z)  + Bu(z) 

- zx(z)  - zx(0)  (24d) 

and 

Z[x*(kT  + T)]  = Gj  Z[x*lkT)]  + Z[u*(kT)]  = GjX*(z)  + u(z) 

J 

= zx*(z)  - zx*(0)  (24e) 

Substituting  G^  and  M-  of  Equations  (18)  and  (19)  in  Eq.  (24e)  and 
!>  b 

expressing  it  in  the  form  of  Eq.  (24d)  yields  ! 

I 

[I  + |2(AT)^]|-f^x*(z)  - (I  - y at  + |y(AT)^]|-^x*(0)  j 

= Ax*(z)  + Bu(z)  {24f)  I 

Comparing  Eq.  (24d)  with  (24f)  yields 

t , 

Z[x(kT)]  s [I  + |y(AT)^]|-g[x(z)  - [I  - 7 at  + •jy(AT)2]|-^x(0)  (24g) 

Note  that  the  approximate  numerical  differentiator  in  Eq.  (23a)  differs  from 
that  of  Eq.  (24g)  with  x(0)  - 0 by  a weighting  factor  [I  + jY(AT)^].  In 
the  same  fashion,  if  G - G^  and  M » are  used,  then  we  have 

Z[x(kT)]  * (I  +^(AT)2]''[I  + •|o^AT)^]^|^x(z) 

-[I  + ^(AT)^]"’[I  - y AT  + l^(AT)^  - T57(AT)^]|-^  x(0)  (24h) 

i 
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In  general,  the  approximate  numerical  differentiator  can  be  written  as 
Z[x(kT)]  s - Dj)"'(N.  + Dj)|-g|x(z)  - AT(Nj  - Dj)’’  Oj  ^^x(0) 

= A(Nj  - Dj)"’(Nj.  + Oj)  f^x(z)  - 2A(Nj  - Dj)"’  ^ x(0)  (241) 

where 

Gj  s Nji  j - 1,3.5,... 

The  approximate  models  obtained  by  using  an  even  number  of  Hj  do  not  have 
ideal  phase  characteristics. 


2.5  SAMPLING  PERIOD 

In  a large  control  system  it  is  often  difficult  to  determine  a 
minimal  common  sampling  period  among  subsystems.  This  method  provides 
more  flexibility  in  determining  the  sampling  period.  In  other  words,  if 
a more  accurate  model  is  used,  a larger  sampling  period  can  be  used. 

Consider  the  commonly  used^  ^ system  matrix  in  Eq.  (22)  and 


substitute  x = T,  or 


Gj  = [I  - j Ax]  * [I  + j Ax] 


(25a) 


To  study  the  relationships  between  the  sampling  period  of  G^  and  that 
of  a more  accurate  model  In  Eq.  (18),  we  equate  G^  with  in  Eq.  (18): 

[I  - j Ax]“’[l  + i-Ax]  » [I  - j AT  + i_(AT)^]‘’[l  + j AT  + ■}y(AT)^] 


(25b) 


With  no  loss  of  generality,  we  can  assume  that  A is  a diagonal  matrix 

with  all  eigenvalues  located  on  the  diagonal.  The  absolute  value 

of  the  largest  eigenvalue,  designated  X , is  used  to  determine  the 

m 

minimal  sampling  period.  One  equation  of  Eq.  (25b)  is 


(1  - ^ x„x)-'(i  * 'j  - n - ^ V * " j V * y'y' 


Solving  Eq.  (25c)  yields 


X - T p =■  < T 

’ TZ^ 


Therefore  a larger  sampling  period  can  be  used  if  instead  of  is 
used  to  approximate  the  system.  In  a similar  manner,  if  G^  were  used. 


, ' * 5o<V''  , , 

X = T 1 IT  < T 


’ * To''J>' 


Likewise,  If  were  used,  then 


X - T 


’ ^ y\p  * 


Note  that  the  sampling  period  x In  Eq.  (25a)  is  always  smaller  than  the 
sampling  period  T In  Equations  (26a),  (26b)  or  (26c). 

2.6  ILLUSTRATIVE  EXAMPLE 

Consider  an  unstable  continuous-time  state  equation: 

x(t)  ■■  Ax(t)  + Bu(t) 

x(0)  - x^  I 

where 


1 2 

A = , B 

3 


, x(0) 


, x(t) 


x^  (t) 
X2(t) 


and  u(t)  Js  the  unit-step  function.  The  approximate  discrete-time 
state  equation  with  sampling  period  T ■ ^ Is  required. 

x*(k  + 1)  » Gx*(k)  + Mu(k) 

Four  approximate  models  are  compared  in  this  example. 

Two  popular  models  often  used  In  Industry  are 

G =.  G*  « I + AT  + ^AT)^ 

M = M*  * [G*  - I]a“’  B - T[|  + j AT]B 


where 


* [’• 


I46875  0.3125 
l»6875  0.6875 


if 

0.625 

0.0625 

and  = 

0.3125 

0.>25 

(27b) 


(28a) 

(28b) 


where 


G = Gg  = I + AT  + '^(AT)^  + -^(AT)^  + '5y(AT)** 


M = Mg  = T[I  + AT  + •^(AT)^  + •^(AT)^]B 


1.456868  0.383138 
0.574707  0.499023 


* r°- 
^ 0. 


648438  O.O5306O 
324219  O.I65039 


(29a) 

(29b) 


Two  other  models,  obtained  from  matrix  continued  fractions  are 


G = Gj  - [I  - j AT]"’  [I  + i AT] 
M = = T[1  - j AT]"’  B 


(30a) 


(30b) 


where 


where 


I.A61538  0.410256  0.666666  0.051282 

and  M_  = 

0.615384  0.435897  ^ 0.333333  0.179487 


G = Gg  - [I  - i AT  + i-(AT)^]"’ [I  + 1 AT  + |y(AT)^] 


M = - T[l  - j AT  + |j(AT)^]"’  B 


1.454246  0.388804  0.648648  0.051968 

G_  = and  = 

^ 0.583206  0.482235  ^ 0.324324  0.168417 


Note  that  the  model  in  Eq.  (30)  is  the  model  used  by  the  referenced 
authors 

The  exact  solution  of  Eq.  (27a)  is 

u\  16  2t  3 -5t  6 

x,(0.y-e  -55-e  -j 


!.■,  8 2t  ^ 3 .-5t  2 

X2(t)  . y . + yy  . - y 


(31b) 


(32b) 


The  responses  at  the  sampling  instants  k ■ 0,1, 2, 3. 4 of  the  exact  model 
and  the  four  approximations  are  shown  In  Table  2-1  [state  x^(kT)]  and 
Table  2-2  [state  X2(kT)].  From  the  tables  we  observe  that  Eq.  (31)  gives 
the  best  approximation. 
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TABLE  2-1 


COMPARISON  OF  STATE  x, (kT) 


EXACT  SOL. 

DIRECT  APPROXIMATION 

MATRIX  CONI 

r.  FRACTION 

Eq.  (32a) 

Eq.  (28) 

Eq.  (29) 

Eq.  (30) 

Eq.  (31) 

1 

1 

1 

1 

1 

2.5'f4 

2.468 

2.541 

2.589 

2.544 

5.006 

4.811 

5.002 

5.145 

5.006 

9.0A2 

8.595 

9.036 

9.380 

9.041 

15.689 

14.731 

15.676 

16.436 

15.686 

TABLE  2-2 

COMPARISON  OF  STATE  X2(kT) 


D 

■ 

EXACT  SOL. 

DIRECT  APPROXIMATION 

MATRIX  CONT.  FRACTION 

Eq.  (32b) 

Eq.  (28) 

Eq.  (29) 

Eq.  (30) 

Eq.  (31) 

D 

0.00 

1 

1 

1 

1 

1 

0.25 

1.558 

1.593 

1.563 

1.564 

1.558 

0.50 

2.728 

2.690 

2.730 

2.788 

2.728 

3 

0.75 

4.728 

4.542 

4.726 

4.894 

4.728 

B 

1.00 

8.046 

7.589 

8.041 

8.419 

8.045 

3.0  STABILITY  OF  COUPLED  MULTIVARIABLE  MISSILE  SYSTEMS 

3.1  INTRODUCTION 


The  accurate  description  of  most  practical  systems,  for  example  both  a 

2 12 
small  semi-active  terminal  homing  missile  system  and  an  aircraft  system  , 

result  in  high  order  coupled  multivariable  differential  equations.  Linear 

representations  of  these  systems  are  by  a set  of  coupled  high-order  differential 

equations  or  a matrix  differential  equation.  A primary  concern  in  the  design 

of  these  multivariable  systems  is  the  stability  problem.  One  conventional 

approach  is  to  formulate  the  system  into  a high  dimensional  state  equation, 

then  to  determine  the  stability  by  either  directly  evaluating  the  roots  of  the 

scalar  characteristics  polynomial,  indirectly  applying  the  Routh  criterion'^, 

11( 

or  application  of  Jury's  inner  theory  on  the  characteristic  polynomial. 
However,  the  determination  of  a characteristic  polynomial  for  a large  dimen- 
sional system  is  tedious.  Moreover,  if  a system  is  modeled  as  a matrix 
differential  equation,  it  is  more  natural  to  determine  the  stability  directly 
from  the  matrix  polynomial  than  indirectly  from  a scalar  polynomial.  Some 
approaches  have  been  proposed  to  determine  the  stability  of  a multivariable 
system  directly  from  the  matrix  polynomial.  Papaconstant inou' ^ suggested  a 
scheme  for  testing  stability  of  polynomial  matrices.  In  his  work,  a 
recursive  algorithm  was  developed  to  compare  the  normalized  largest  eigen- 
values with  unity.  However,  the  method  requires  the  calculation  of  the 
eigenvalues  of  largest  moduli.  This  procedure  is  difficult  due  to  the  problem 
of  convergence  of  the  eigenvalue  algorithm.  Recently,  Shleh  and  Sacheti'^ 


partially  extended  the  scalar  Routh  criterion  ^ to  the  matrix  case.  In  this 
work,  It  is  shown  that,  If  a matrix  polynomial  B(s)  * Is*^  + ' + ...  + Bj 

is  given,  a matrix  Routh  array  can  be  constructed  by  using  the  following 
matrix  Routh  algorithm: 


'•j'Wzj  « 


where  I 


■=■  + 1 n even 


n odd 


^2,j  “ ®n+2-2j  j ■ 


~ ''i-2,j+1  ■ ^l-2^i-l,j+l  ‘ “ 1.2,...,  j = 3,^,... 

”l  “ *"1 ,1  ^^i+1 ,1^  * “ 1,2,...,  n 

det  (C.^,  j)  0 (1) 

A sufficient  condition  for  stability  of  the  det  [B(s)]  is  that  all  the  "matrix 
quotients"  Hj  be  real,  symmetric,  positive  definite  matrices.  Note  that  this 
sufficient  condition  deals  only  with  H,  and  not  C.  , (the  block  elements  in 
the  first  column  of  the  matrix  Routh  array).  Liapunov  theory  with  the  state 
equation  in  the  controllable  block  companion  (controllable  phase-variable) 
form  was  used  to  derive  their  result. 

In  this  report,  we  develop  two  approaches  for  determining  the  stability  of 
a class  of  multivariable  systems.  One  approach  uses  the  "matrix  quotients" 
that  are  developed  from  an  alternate  matrix  Routh  algorithm  and  a state 
equation  in  the  observable  block  companion  form'^.  The  other  approach  uses 
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the  block  elements  In  the  first  column  of  the  matrix  Routh  array.  Two 
sufficient  conditions  and  three  necessary  conditions  are  derived  for  the 
stability  of  matrix  polynomials,  thereby  partially  extending  the  scalar 
Routh  criterion  to  the  matrix  Routh  criterion. 

3.2  SUFFICIENT  CONDITIONS 

The  objective  of  this  report  is  to  establish  the  criteria  for  the  stability 
of  the  following  matrix  state  equations. 

i-i 

Y B.D  x(t)  = [0],  B , = I (2a) 

1-1  ' " ' 


and 

D*  'x(O)  - • “ 1,2, 3,.. .,n  (2b) 

where  x(t)  Is  the  m-d Imens ional  state  vector.  Bj,  I,  and  [O]  are  m x m 
real  constant  matrix,  identity  matrix  and  null  matrix  respectively.  For  the 
scalar  case,  it  Is  well  known  that  a system  is  asymptotically  stable  if  and 

only  if  the  Routh  array  elements  in  the  first  column  are  all  positive.  Shieh 

16  13 

and  Sacheti  partially  extended  the  Routh  criteria  to  the  matrix  case  and 

derived  a sufficient  condition  for  the  stability  of  a multivariable  system 

in  Eq.  (2)  from  the  controllable  block  companion  form.  In  this  report  we 

derive  some  sufficient  and  some  necessary  conditions  for  the  system  in  Eq. 

(2)  from  the  observable  block  companion  form. 
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Let  us  rewrite  the  system  In  Eq.  (2)  Into  the  following  observable  block 


companion  form; 


[k]  = [B][x] 
U(0)]  » [a] 


where 


0 0 0 . 0 -B, 


I 0 0 . 0 -B^ 


[B]  = 0 I 0 . 0 -B, 


0 0 0.  I -B 


The  dimensions  of  the  matrix  [B] , the  block  elements  Bj,  and  state  vector  [x] 
are  (nxm)  x (nxm) , m x m,  and  (nxm)  x 1 respectively.  Equation  (3)  can  be 
transformed  into  the  block  Schwarz  form  by  using  the  following  linear  trans- 


formation: 


[x]  - [K,]Iy] 


[y]  = [K,l‘’[B][K,]ty]  - [Ally] 


7 


where 


[K,] 


• 1 

-1 

-1 

- 

1 . 

D D ' 

Vl,2Vl, 

• V3.3n-3.l 

0 

° 

.1 

-1 

.1 

0 . 

1 

0 

-1 

0 

0 . 

0 

0 

■’52°5i 

0 

-1 

'’33°3l 

0 

0 . 

0 

1 

r 

0 

0 

“23“2l 

0 . 

0 

0 1 

1 

1 

1 

0 

r 

“32°3l 

0 

0 . 

0 

0 j 

1 

0 

• 1 

1 

1 

0 

'>22“2i 

0 . 

0 

0 j 

1 

0 

1 ^ 

( 

1 1 

1 

1 

0 

1 

0 . 

0 

0 j 

0 

1 ° 

1 ° 

! ' J 

and 


[A]  = 


0 -A, 


0 

I 

0 

0 


0 

-A. 


0 

0 


0 
0 
0 

0 -A 

I 

1 -A 


0 

0 

0 

h-1 


(Ac) 


(Ad) 


The  dimension  of  each  block  element  In  [A]  and  [K^]  Is  m x m.  The  block 
elements  Dj  j,  having  dimension  m x m,  In  Eq.  (Ac)  can  be  obtained  from  the 
following  alternate  matrix  Routh  algorithm  and  alternate  matrix  Routh  array 
which  are  different  from  those  In  Eq.  (1). 


where 


“ *^1-2, j+1  ■ '’i-l,j+l^i-2  j “ ‘ ■ 3.^.-.- 

" 1.2, ...,n 

det  [0,^,  ,]  ^ 0 (5c) 

Shleh  and  Sachetl'^  have  shown  that  if  H.  = C.  ^ for  i = 1,2,...,n 

In  Eq.  (1)  are  positive  definite,  then  the  system  In  Eq.  (2)  Is  asymptotically 
stable.  Here,  we  show  similar  results  when  replacing  H.  by  M.. 

Theorem  I ; 

If  {M.}  i = 1,2,...,n  in  Eq.  (5)  are  positive  definite,  then  the  system 
In  Eq.  (2)  Is  asymptotically  stable. 


Proof: 


Performing  the  following  new  transformation 


on  Eq.  (4)  yields 


[y]  = [K2][z] 


[i]  = [K2]''[A][K2][z] 


= [F][z] 


where 


0 0 


0 0 


n-1 ,1 


0 0 


D.,  0 


0 0, 


0 -M 


m"’,  0 -m"\ 

n-1  n-1 


0 M 


[F]  = . 
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The  Ifnear  transformation  matrix  Ik]  between  x coordinates  and  z coordin- 


ates Is 


W = [Kllz]  - IK,][K2][z] 


Now,  consider  the  following  quadratic  equation: 

V - [q]^[P][q] 


where 


M 0 
n 


0 0 


0 M 


0 0 


0 0 


M^  0 


0 0 


0 M, 


and  T in  Eq.  (8a)  designates  transpose. 
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since  {M.}  are  positive  definite  which  implies  that  P is  positive  definite, 
V Is  positive  definite.  The  derivative  of  V is 

V = [q]''^[PF  + F''’p][q] 

= -[q]^[Q][q]  Oa) 

where 


0 0 


0 0 


0 0 


[Cl]  = . 


0 0 


0 0 


0 0 


0 0 


0 21 


From  Equations  (8)  and  (9)  we  can  see  that  V is  a Liapunov  function.  Hence, 
we  conclude  that  the  system  in  Eq.  (2)  Is  asymptotically  stable. 

From  the  result  obtained  in  Theorem  I,  we  establish  another  sufficient 
condition  for  the  stability  of  the  system  in  Eq.  (2)  by  using  the  block 
elements  Dj  ^ In  the  matrix  Routh  array  in  Eq.  (5)  Instead  of  the  M.  in 
Eq.  (5). 

Theorem  2; 

If  {Oj  1 = are  positive  definite,  the  eigenvalues  of 

{D.  1 = 1,3,5 are  positive  and  real,  and  {Dj  i = 1,3,5 

^*^1+1  1*^1  “ 2*^.6,...,  are  Hermitlan,  the  system  [Eq.  (2)]  is  asymptotically 

stable. 

In  order  to  prove  Theorem  2,  we  need  the  following  lemma  which  is  due  to 
KyFan’®  [P.  137]. 


t 


i 

Lemma  1 . Let  Kj  be  positive  definite  and  K2  such  that  ^1^2  Hermltian.  j 

Then  KjK2  is  positive  definite  if  and  only  if  the  eigenvalues  of  K2  are 
positive  and  real.  In  the  following  lemma,  we  switch  the  conditions  on 
Kj  and  K2  yielding  the  same  result. 

Lemma  2.'  Let  K2  be  positive  definite  and  Kj  such  that  KjK2  is  Hermltian. 

i 

The  K^K2  positive  definite  if  and  only  if  the  eigenvalues  of  Kj  are  i 

positive  and  real . 

Proof;  Since  K2  is  positive  definite  which  implies  is  positive 
definite,  where  T designates  transpose,  it  is  seen  from  lemma  1 that 
is  positive  definite  if  and  only  if  the  eigenvalues  of  k|  are 
positive  and  real.  But  = (K^K2)^;  i.e.,  K^K2  is  positive  definite 

if  and  only  if  the  eigenvalues  of  Kj  are  positive  and  real. 

Lemma  3.  If  K2  is  positive  definite  and  K^K2  is  symmetric,  then 
is  symmetric. 

Proof:  Since  {KjK2)^  = ^2*^1  “ * *^1*^2  implies  k|  = K2'KjK2. 

Hence  (k”’k2)^  = K2(Kj’)^  = K2K2’k“’k2  = k'’k2;  i.e.,  Kj’K2  is  symmetric. 

Proof  of  Theorem  2; 

By  lemma  3t  we  know  that  jD.  ^ Is  symmetric  for  i = 1,2,...  . By 

lemma  2 or  3,  we  know  that  Mj  = ^Oj  ^ is  positive  definite.  Hence,  the 

system  in  Eq.  (2)  is  asymptotically  stable  following  the  results  of  Theorem  1. 

In  order  to  show  an  application  of  Theorem  1 and  Theorem  2,  let  us  consider 
the  following  matrix  characteristic  equation: 

As^  + Bs  + C » 0 (10a)  1 


! 

t. 
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where 


I 0 

3 I 

, B = 

and  C = 

0 1 

1 3 

\ 

o 

N> 

1 

Routh  algorithm  of  Eq.  (1),  we  obtain 
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3 
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In  this  case,  no  conclusion  can  be  drawn  from  the  sufficient  condition 
established  by  Shleh  and  Sachetl'^.  However,  If  we  arranged  the  matrices  A, 
B and  C according  to  Eq.  (5),  we  have 


From  Theorem  1,  we  see  that  the  system  is  asymptotically  stable. 

This  example  shows  the  application  of  Theorem  2.  Let  us  consider  the 
following  matrix  characteristic  equation: 

+ B,S  + = 0 (Ila) 

where 


In  Bel Iman’®  [p.67,  p.lOl]  it  is  shown  that,  if  A,,  B^  and  are  positive 

2 

definite,  then  the  roots  of  det  [AjS  + BjS  + Cj]  •»  0 have  negative  real 
parts.  But  in  this  example,  no  conclusion  can  be  made  from  Bellman's  results. 


,1 


However,  we  know  that 


and 


3 1.63 

1.63  0.9 
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D 


31 


(11b) 


which  are  synmetrlc,  Is  positive  definite,  and  the  eigenvalues  of  Aj  and 
Cj  are  positive  and  real.  Therefore,  from  Theorem  2 we  conclude  that  the 
system  In  Eq.  (I  la)  is  asymptotically  stable. 

3.3  NECESSARY  CONDITIONS 

In  this  section  we  establish  some  necessary  conditions  for  the  stability 
of  multivariable  systems.  The  failure  to  satsify  the  necessary  conditions 
for  stability  Is  equivalent  to  the  sufficient  conditions  for  the  Instability 
of  the  same  systems;  i.e.. 

Theorem  3. 

If  {Mj}  i = 1,2 n are  symmetric  such  that,  there  exists  one  {M.} 

I = I,2,...,n  which  is  negative  definite,  negative  semi-definite,  or 
indefinite,  then  the  system  in  Eq.  (2)  is  unstable. 

Proof; 

Suppose  the  system  is  asymptotically  stable  and  one  of  {Mj}  is  negative 
definite,  negative  semi-definite,  or  indefinite.  Since  the  stability  is 
invariant  under  the  linear  transformation  and  the  matrix  F in  Eq.  (7)  is 
a stable  matrix.  Let  us  consider  the  following  equation: 

XF  + F^X  = -Q  (12) 

1 Q 

where  Q Is  a matrix  defined  in  Eq.  (9b).  By  Theorem  k in  Bellman  [p.  239] 
we  know  that  Eq.  (12)  has  a unique  solution.  Since  Q is  positive  semi-definite 
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we  conclude  that  the  solution  X of  Eq.  (12)  is  also  positive  semi -definite. 

It  is  easy  to  verify  that  the  matrix  P which  was  defined  In  Eq.  (8b) 
satisfies  Eq.  (12).  Therefore  X = P is  positive  semi -definite.  This  implies 


that  at  least  one  of  the  {Mj}  Is  positive  seml-def Ini te  and  others  positive 
definite.  This  contradicts  our  assumption  that  one  of  the  {Mj}  is  negative 
definite,  negative  seml-def Inite,  or  Indefinite.  Hence  the  system  In  Eq.  (2) 
Is  unstable  If  one  of  the  {Mj}  Is  negative  definite,  negative  seml-def inite. 


or  indefinite. 


To  show  an  application  of  Theorem  3,  consider  the  example 
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Is  symmetric  and  indefinite.  According  to  Theorem  3,  the  system  is 
unstable. 
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The  following  theorem  is  another  criteria  for  an  unstable  system. 


Theorem  k. 

If  = I and  the  trace  of  ‘'21  = «n  is  negative,  then  the 

system  in  Eq.  (2)  Is  unstable,  where  and  D2^  are  defined  in  Eq.  (5a). 
Proof; 

The  matrix  characteristic  equation  of  the  system  in  Eq.  (2)  is 
[D(S)]  - OjiS"  * Dj.s"-'  ... 

where  1=1  if  n is  even  and  i =2  if  n Is  odd. 

Since  the  sum  of  the  eigenvalues  of  the  system  in  Eq.  (2)  is  equal  to  the 
negative  value  of  the  trace  of  , this  implies  that  there  exists  some 
eigenvalues  of  det  [D(S)]  which  are  positive.  Hence,  the  system  is  unstable. 

The  next  criteria  Is  another  necessary  condition  which  we  state  as 
fol lows. 

Theorem  5. 

If  det  D,,  > 0 and  det  D.  < 0,  or  det  D,,  < 0 and  det  D,  > 0,  and 
11  I ,n  II  i ,n  ’ 

Dll,  Dj  are  defined  In  Eq.  (l^i)  then  the  system  in  Eq.  (2)  is  unstable. 
Proof; 

Since  the  system  in  Eq.  (2)  has  the  matrix  characteristic  equation 
[D(S)]  In  Eq.  (14),  then  we  expand  the  det  [D(S)].  We  find  the  constant 
term  is  equal  to  det  B,  ■ det  D,  If  det  D,,  > 0 and  det  D.  „ < 0,  this 
implies  that  the  coefficient  of  the  polynomial  det  [D(S)]  has  a negative  sign. 
We  can  then  conclude  that  the  det  [D(S)]  = 0 
real  part.  Hence  the  system  Is  unstable. 


4.0  CONCLUSIONS 


A method,  based  on  a mode)  reduction  technique,  for  accurately 
representing  continuous-time  state  equations  by  discrete-time  state 
equations  has  been  presented.  Matrix  continued  fractions  have  been  used 
as  a basis  for  the  derivation.  The  relationships  among  commonly  used 

1 

i approximation  models  and  the  proposed  models  have  been  discussed.  It 

j has  been  shown  that  the  models  developed  by  Tustin,  Boxer  and  Thaler  and 

j Shieh  et.al.  are  special  cases  of  the  methods  presented  here.  Also,  it 

I has  been  shown  that  if  a more  accunjte  model  is  used,  a larger  sampling 

I period  can  be  applied.  Thus  the  o'fficulty  of  selecting  a common 

I 

I 

sampling  period  among  the  subsystems  of  the  T6  missile  can  be  reduced. 

Some  necessary  and  sufficient  conditions  have  been  developed  for  the 
stability  of  multivariable  systems.  A linear  block  transformation  has 
been  derived  for  transforming  the  coordinates  of  an  observable  block 
companion  form  to  the  coordintes  of  a block  Schwart  form.  The  proposed 
method  has  partially  extended  the  scalar  Routh  criterion  to  the  matrix 
Routh  criterion  to  a class  of  multivariable  systems. 

Thus  each  subsystem  of  the  T6  missile  can  be  more  accurately  repre- 
sented by  the  discrete-time  models  proposed  in  this  report.  Also,  the 
stability  of  the  coupled  yaw  and  roll  system  can  be  determined  by  the 
proposed  approach.  The  proposed  methods  in  this  report  can  be  applied 
in  general  to  missile  systems  of  the  T6  missile  class. 
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